###Simulating Trajectories
tmax <- 10
X1_0 <- 600
X2_0 <- 30
X3_0 <-10^5
rho <- 0.108
delta <- 0.5
eta <- 9.5*10^(-6)
lambda <- 36
N1 <- 1000
C <- 3
sigma1 <- 0.1
sigma2 <- 0.1
sigma2 <- 0.1
make_path <- function(tmax, x1_0, x2_0, x3_0, rho, delta, eta, lambda, N1, C, sigma1, sigma2, sigma3, deltat, n){
#Initializing
nosteps <- tmax/deltat
x1 <- numeric(nosteps + 1)
x2 <- numeric(nosteps + 1)
x3 <- numeric(nosteps + 1)
x <- matrix(data = c(x1, x2, x3), ncol = 3)
x[1,] <- c(x1_0, x2_0, x3_0)
sigma <- c(sigma1, sigma2, sigma3)
#Simulationg the driving BM
dW1 <- rnorm(n = nosteps, mean = 0, sd = sqrt(deltat))
dW2 <- rnorm(n = nosteps, mean = 0, sd = sqrt(deltat))
dW3 <- rnorm(n = nosteps, mean = 0, sd = sqrt(deltat))
dW <- matrix(data = c(dW1, dW2, dW3), ncol = 3)
#Defining drift function
drift <- function(x){
x1 <- x[1]
x2 <- x[2]
x3 <- x[3]
y1 <- lambda - rho*x1 - eta*x1*x3
y2 <- eta*x1*x3 - delta*x2
y3 <- N1*delta*x2 - C*x3
y <- c(y1, y2, y3)
y
}
#Iterating
for (k in (1:(nosteps))){
x[k+1,] <- x[k,] + drift(x[k,])*deltat + sigma*x[k,]*dW[k,] + 0.5*sigma*x[k,]*sigma*(dW[k,]*dW[k,] - deltat)
}
#Creating output and adding time variable
t <- (0:nosteps)*deltat
out <- data.frame(x, t)
colnames(out) <- c("x1", "x2", "x3", 't')
out
}
subsample <- function(x){
x <- as.matrix(x)
x_fine <- matrix(data = 0, nrow = 500, ncol = 4)
for (i in (1:500)){
x_fine[i,] <- x[20*(i-1) + 1,]
}
x_fine <- as.data.frame(x_fine)
colnames(x_fine) <- c('x1','x2','x3','t')
x_coarse <- matrix(data = 0, nrow = 100, ncol = 4)
for (i in (1:100)){
x_coarse[i,] <- x[100*(i-1) + 1,]
}
x_coarse <- as.data.frame(x_coarse)
colnames(x_coarse) <- c('x1','x2','x3','t')
out <- list(x_fine, x_coarse)
out
}
path <- make_path(tmax = 10, x1_0 = 600, x2_0 = 30, x3_0 = 10^5, rho = 0.108, delta = 0.5, eta = 9.5*10^(-6), lambda = 36, N1 = 1000, C = 3, sigma1 = 0.1, sigma2 = 0.1, sigma3 = 0.1, deltat = 0.001)
sampled_path <- subsample(path)
x_fine <- sampled_path[[1]]
x_coarse <- sampled_path[[2]]
###Plotting the trajectory
path_plot_fine <- plot_ly(x_fine, x = ~x1, y = ~x2, z = ~x3,
type = 'scatter3d', mode = 'lines+markers',
line = list(width = 6, color = ~t, colorscale = 'Viridis'),
marker = list(size = 3.5, color = ~t, colorscale = 'Greens', cmin = -20, cmax = 50))
path_plot_fine
path_plot_coarse <- plot_ly(x_coarse, x = ~x1, y = ~x2, z = ~x3,
type = 'scatter3d', mode = 'lines+markers',
line = list(width = 6, color = ~t, colorscale = 'Viridis'),
marker = list(size = 3.5, color = ~t, colorscale = 'Greens', cmin = -20, cmax = 50))
path_plot_coarse
Euler-Maruyama Estimator
For making inference, we need to have a likelihood for a ny given
parameter \(\theta\). This will be the
sum of the likelihoods of each step, under the EM-hypothesis. Thus, for
each increment, we need to assign a likelihood of that increment, given
the parameter vector. In this case, the stochastic increments to each
coordinate are independent, due to the diagonal structure of \(\Sigma\). So the stochastic quantities here
are Normal and Chi-square distributed. But we have a sum here, but one
is a deterministic transformation of the other. But there might be more
that one dW that give the same increment. In that case i suppose we have
to assign both the likelihoods. Do we need to take a convolution to find
the density of the sum? Apparently, we are to use a Gaussian
approximation as given on the slides. The books probably cover this as
well. Upon looking at the slides again, the EM-based estimator seems
quite doable. Implementing likelihood evaluation for single step
likelihood_step <- function(theta, x, x_next, deltat){
#Defining tunable and constant parameters
lambda <- theta[1]
N1 <- theta[2]
C <- theta[3]
sigma1 <- theta[4]
sigma2 <- theta[5]
sigma3 <- theta[6]
#Defining drift function
drift <- function(x){
x1 <- x[[1]]
x2 <- x[[2]]
x3 <- x[[3]]
y <- numeric(length = 3)
y[1] <- lambda - rho*x1 - eta*x1*x3
y[2] <- eta*x1*x3 - delta*x2
y[3] <- N1*delta*x2 - C*x3
y
}
# need to assign subparameters of theta to the names used in the code
z <- numeric(3)
z[1] <- x[[1]]
z[2] <- x[[2]]
z[3] <- x[[3]]
#z_next <- numeric(3)
#z_next[1] <- x_next
#browser()
y <- drift(x)
a <- log(((2*pi*deltat)^3) * (z[1]*z[2]*z[3]*sigma1*sigma2*sigma3)^2)
b <- ((x_next - x + y*deltat)^2)*c((z[1]*sigma1)^(-2), (z[2]*sigma2)^(-2), (z[3]*sigma3)^(-2))/deltat
c <- b[[1]]*b[[2]]*b[[3]]
a+c
}
#Test chunk
theta <- c(lambda, N1, C, sigma1, sigma2, sigma3)
likelihood_step(theta = theta, x = x_coarse[1,], x_next = x_coarse[2,], deltat = 0.1)
[1] 119262182
#We get a number, so that's good.
Implementing likelihood for entire process
likelihood_total <- function(theta, process, deltat){
nosteps <- length(process[,1])-1
cumulative_likelihood <- 0
for (i in (1:nosteps)){
cumulative_likelihood <- cumulative_likelihood + likelihood_step(theta = theta, x = process[i,], x_next = process[i+1,], deltat = deltat)
}
cumulative_likelihood
}
Testing total likelihood function
likelihood_total(theta = theta, process = x_coarse, deltat = 0.1)
[1] 125389127
#We get a number, which is promising
likelihood_total(theta = theta_0, process = x_coarse, deltat = 0.1)
[1] 5108.032
We see that evaluating the total likelihood function returns the same
number for two different set of parameters. This indicates that that the
total likelihood fucntion doesn’t in fact use the parameter provided for
evaluation. It seems that likelihood_step takes its parameters from the
global envieronment rather than its arguement. Ths would explain why the
likelihood is evaluated to the same for any parameter. Now we got
different numbers.
Now, we try to optimize this for some parameter vector theta.
theta_0 <- c(50, 2000, 10, 0.5, 0.5, 0.5)
optim(par = theta_0, fn = likelihood_total, process = x_coarse, deltat = 0.1, lower = c(0,0,0,0,0,0), method = "L-BFGS-B")
Error in optim(par = theta_0, fn = likelihood_total, process = x_coarse, :
L-BFGS-B needs finite values of 'fn'
It would appear that i have made a mistake: Optim outputs the initial
parameters, and states that convergence was not reached. Another mistake
is not finishing a coding project before taking a 2-week break: It’s a
little difficult to jumo in again.
It takes some time to run now, and we don’t get convergence or the
right results. The total likelihood function was called over 400 times,
probably because of lack of convergence, which might also be why it
takes so long. Question: are we maximizing or minimizing? What is the
parameter space? What exactly is the quantity being computed in the step
likelihood function? I think it is the likelihood, formula from the
slides. Let’s try with some parameters that are closer to the real ones.
We get completely atrocious estimates. We try bounding barameter values
from below with zero and using L-BFGS-B method. It seems that we get
infinite values from total likelihood, which is unfortunate. Maybe the
log is fucking me?
LS0tDQp0aXRsZTogIklTREUgRXhlcmNpc2VzIFdlZWsgNCINCm91dHB1dDogaHRtbF9ub3RlYm9vaw0KLS0tDQoNCiMjI1NpbXVsYXRpbmcgVHJhamVjdG9yaWVzDQoNCmBgYHtyfQ0KdG1heCA8LSAxMA0KWDFfMCA8LSA2MDANClgyXzAgPC0gMzANClgzXzAgPC0xMF41DQpyaG8gPC0gMC4xMDgNCmRlbHRhIDwtIDAuNQ0KZXRhIDwtIDkuNSoxMF4oLTYpDQpsYW1iZGEgPC0gMzYNCk4xIDwtIDEwMDANCkMgPC0gMw0Kc2lnbWExIDwtIDAuMQ0Kc2lnbWEyIDwtIDAuMQ0Kc2lnbWEzIDwtIDAuMQ0KYGBgDQoNCmBgYHtyfQ0KbWFrZV9wYXRoIDwtIGZ1bmN0aW9uKHRtYXgsIHgxXzAsIHgyXzAsIHgzXzAsIHJobywgZGVsdGEsIGV0YSwgbGFtYmRhLCBOMSwgQywgc2lnbWExLCBzaWdtYTIsIHNpZ21hMywgZGVsdGF0LCBuKXsNCiAgI0luaXRpYWxpemluZw0KICBub3N0ZXBzIDwtIHRtYXgvZGVsdGF0DQogIHgxIDwtIG51bWVyaWMobm9zdGVwcyArIDEpDQogIHgyIDwtIG51bWVyaWMobm9zdGVwcyArIDEpDQogIHgzIDwtIG51bWVyaWMobm9zdGVwcyArIDEpDQogIHggPC0gbWF0cml4KGRhdGEgPSBjKHgxLCB4MiwgeDMpLCBuY29sID0gMykNCiAgeFsxLF0gPC0gYyh4MV8wLCB4Ml8wLCB4M18wKQ0KICBzaWdtYSA8LSBjKHNpZ21hMSwgc2lnbWEyLCBzaWdtYTMpDQogIA0KICAjU2ltdWxhdGlvbmcgdGhlIGRyaXZpbmcgQk0NCiAgZFcxIDwtIHJub3JtKG4gPSBub3N0ZXBzLCBtZWFuID0gMCwgc2QgPSBzcXJ0KGRlbHRhdCkpDQogIGRXMiA8LSBybm9ybShuID0gbm9zdGVwcywgbWVhbiA9IDAsIHNkID0gc3FydChkZWx0YXQpKQ0KICBkVzMgPC0gcm5vcm0obiA9IG5vc3RlcHMsIG1lYW4gPSAwLCBzZCA9IHNxcnQoZGVsdGF0KSkNCiAgZFcgPC0gbWF0cml4KGRhdGEgPSBjKGRXMSwgZFcyLCBkVzMpLCBuY29sID0gMykNCiAgDQogICNEZWZpbmluZyBkcmlmdCBmdW5jdGlvbg0KICBkcmlmdCA8LSBmdW5jdGlvbih4KXsNCiAgICB4MSA8LSB4WzFdDQogICAgeDIgPC0geFsyXQ0KICAgIHgzIDwtIHhbM10NCiAgICB5MSA8LSBsYW1iZGEgLSByaG8qeDEgLSBldGEqeDEqeDMNCiAgICB5MiA8LSBldGEqeDEqeDMgLSBkZWx0YSp4Mg0KICAgIHkzIDwtIE4xKmRlbHRhKngyIC0gQyp4Mw0KICAgIHkgPC0gYyh5MSwgeTIsIHkzKQ0KICAgIHkNCiAgfQ0KICANCiAgI0l0ZXJhdGluZw0KICBmb3IgKGsgaW4gKDE6KG5vc3RlcHMpKSl7DQogICAgeFtrKzEsXSA8LSB4W2ssXSArIGRyaWZ0KHhbayxdKSpkZWx0YXQgKyBzaWdtYSp4W2ssXSpkV1trLF0gKyAwLjUqc2lnbWEqeFtrLF0qc2lnbWEqKGRXW2ssXSpkV1trLF0gLSBkZWx0YXQpDQogIH0NCiAgDQogICNDcmVhdGluZyBvdXRwdXQgYW5kIGFkZGluZyB0aW1lIHZhcmlhYmxlDQogIHQgPC0gKDA6bm9zdGVwcykqZGVsdGF0DQogIG91dCA8LSBkYXRhLmZyYW1lKHgsIHQpDQogIGNvbG5hbWVzKG91dCkgPC0gYygieDEiLCAieDIiLCAieDMiLCAndCcpDQogIG91dA0KfQ0KYGBgDQoNCmBgYHtyfQ0Kc3Vic2FtcGxlIDwtIGZ1bmN0aW9uKHgpew0KICB4IDwtIGFzLm1hdHJpeCh4KQ0KICANCiAgeF9maW5lIDwtIG1hdHJpeChkYXRhID0gMCwgbnJvdyA9IDUwMCwgbmNvbCA9IDQpDQogIGZvciAoaSBpbiAoMTo1MDApKXsNCiAgICB4X2ZpbmVbaSxdIDwtIHhbMjAqKGktMSkgKyAxLF0NCiAgfQ0KICB4X2ZpbmUgPC0gYXMuZGF0YS5mcmFtZSh4X2ZpbmUpDQogIGNvbG5hbWVzKHhfZmluZSkgPC0gYygneDEnLCd4MicsJ3gzJywndCcpDQogIA0KICB4X2NvYXJzZSA8LSBtYXRyaXgoZGF0YSA9IDAsIG5yb3cgPSAxMDAsIG5jb2wgPSA0KQ0KICBmb3IgKGkgaW4gKDE6MTAwKSl7DQogICAgeF9jb2Fyc2VbaSxdIDwtIHhbMTAwKihpLTEpICsgMSxdDQogIH0NCiAgeF9jb2Fyc2UgPC0gYXMuZGF0YS5mcmFtZSh4X2NvYXJzZSkNCiAgY29sbmFtZXMoeF9jb2Fyc2UpIDwtIGMoJ3gxJywneDInLCd4MycsJ3QnKQ0KICANCiAgb3V0IDwtIGxpc3QoeF9maW5lLCB4X2NvYXJzZSkNCiAgb3V0DQp9DQpgYGANCmBgYHtyfQ0KcGF0aCA8LSBtYWtlX3BhdGgodG1heCA9IDEwLCB4MV8wID0gNjAwLCB4Ml8wID0gMzAsIHgzXzAgPSAxMF41LCByaG8gPSAwLjEwOCwgZGVsdGEgPSAwLjUsIGV0YSA9IDkuNSoxMF4oLTYpLCBsYW1iZGEgPSAzNiwgTjEgPSAxMDAwLCBDID0gMywgc2lnbWExID0gMC4xLCBzaWdtYTIgPSAwLjEsIHNpZ21hMyA9IDAuMSwgZGVsdGF0ID0gMC4wMDEpDQpzYW1wbGVkX3BhdGggPC0gc3Vic2FtcGxlKHBhdGgpDQp4X2ZpbmUgPC0gc2FtcGxlZF9wYXRoW1sxXV0NCnhfY29hcnNlIDwtIHNhbXBsZWRfcGF0aFtbMl1dDQpgYGANCg0KIyMjUGxvdHRpbmcgdGhlIHRyYWplY3RvcnkNCmBgYHtyIGluY2x1ZGU9RkFMU0V9DQpsaWJyYXJ5KHBsb3RseSkNCmBgYA0KYGBge3J9DQpwYXRoX3Bsb3RfZmluZSA8LSBwbG90X2x5KHhfZmluZSwgeCA9IH54MSwgeSA9IH54MiwgeiA9IH54MywNCiAgICAgICAgICAgICAgICAgICAgICAgICAgdHlwZSA9ICdzY2F0dGVyM2QnLCBtb2RlID0gJ2xpbmVzK21hcmtlcnMnLA0KICAgICAgICAgICAgICAgICAgICAgICAgICBsaW5lID0gbGlzdCh3aWR0aCA9IDYsIGNvbG9yID0gfnQsIGNvbG9yc2NhbGUgPSAnVmlyaWRpcycpLA0KICAgICAgICAgICAgICAgICAgICAgICAgICBtYXJrZXIgPSBsaXN0KHNpemUgPSAzLjUsIGNvbG9yID0gfnQsIGNvbG9yc2NhbGUgPSAnR3JlZW5zJywgY21pbiA9IC0yMCwgY21heCA9IDUwKSkNCnBhdGhfcGxvdF9maW5lDQpgYGANCmBgYHtyfQ0KcGF0aF9wbG90X2NvYXJzZSA8LSBwbG90X2x5KHhfY29hcnNlLCB4ID0gfngxLCB5ID0gfngyLCB6ID0gfngzLA0KICAgICAgICAgICAgICAgICAgICAgICAgICB0eXBlID0gJ3NjYXR0ZXIzZCcsIG1vZGUgPSAnbGluZXMrbWFya2VycycsDQogICAgICAgICAgICAgICAgICAgICAgICAgIGxpbmUgPSBsaXN0KHdpZHRoID0gNiwgY29sb3IgPSB+dCwgY29sb3JzY2FsZSA9ICdWaXJpZGlzJyksDQogICAgICAgICAgICAgICAgICAgICAgICAgIG1hcmtlciA9IGxpc3Qoc2l6ZSA9IDMuNSwgY29sb3IgPSB+dCwgY29sb3JzY2FsZSA9ICdHcmVlbnMnLCBjbWluID0gLTIwLCBjbWF4ID0gNTApKQ0KcGF0aF9wbG90X2NvYXJzZQ0KYGBgDQoNCiMjIyBFdWxlci1NYXJ1eWFtYSBFc3RpbWF0b3INCkZvciBtYWtpbmcgaW5mZXJlbmNlLCB3ZSBuZWVkIHRvIGhhdmUgYSBsaWtlbGlob29kIGZvciBhIG55IGdpdmVuIHBhcmFtZXRlciAkXHRoZXRhJC4gVGhpcyB3aWxsIGJlIHRoZSBzdW0gb2YgdGhlIGxpa2VsaWhvb2RzIG9mIGVhY2ggc3RlcCwgdW5kZXIgdGhlIEVNLWh5cG90aGVzaXMuIFRodXMsIGZvciBlYWNoIGluY3JlbWVudCwgd2UgbmVlZCB0byBhc3NpZ24gYSBsaWtlbGlob29kIG9mIHRoYXQgaW5jcmVtZW50LCBnaXZlbiB0aGUgcGFyYW1ldGVyIHZlY3Rvci4gSW4gdGhpcyBjYXNlLCB0aGUgc3RvY2hhc3RpYyBpbmNyZW1lbnRzIHRvIGVhY2ggY29vcmRpbmF0ZSBhcmUgaW5kZXBlbmRlbnQsIGR1ZSB0byB0aGUgZGlhZ29uYWwgc3RydWN0dXJlIG9mICRcU2lnbWEkLiBTbyB0aGUgc3RvY2hhc3RpYyBxdWFudGl0aWVzIGhlcmUgYXJlIE5vcm1hbCBhbmQgQ2hpLXNxdWFyZSBkaXN0cmlidXRlZC4gQnV0IHdlIGhhdmUgYSBzdW0gaGVyZSwgYnV0IG9uZSBpcyBhIGRldGVybWluaXN0aWMgdHJhbnNmb3JtYXRpb24gb2YgdGhlIG90aGVyLiBCdXQgdGhlcmUgbWlnaHQgYmUgbW9yZSB0aGF0IG9uZSBkVyB0aGF0IGdpdmUgdGhlIHNhbWUgaW5jcmVtZW50LiBJbiB0aGF0IGNhc2UgaSBzdXBwb3NlIHdlIGhhdmUgdG8gYXNzaWduIGJvdGggdGhlIGxpa2VsaWhvb2RzLiBEbyB3ZSBuZWVkIHRvIHRha2UgYSBjb252b2x1dGlvbiB0byBmaW5kIHRoZSBkZW5zaXR5IG9mIHRoZSBzdW0/IEFwcGFyZW50bHksIHdlIGFyZSB0byB1c2UgYSBHYXVzc2lhbiBhcHByb3hpbWF0aW9uIGFzIGdpdmVuIG9uIHRoZSBzbGlkZXMuIFRoZSBib29rcyBwcm9iYWJseSBjb3ZlciB0aGlzIGFzIHdlbGwuDQpVcG9uIGxvb2tpbmcgYXQgdGhlIHNsaWRlcyBhZ2FpbiwgdGhlIEVNLWJhc2VkIGVzdGltYXRvciBzZWVtcyBxdWl0ZSBkb2FibGUuDQpJbXBsZW1lbnRpbmcgbGlrZWxpaG9vZCBldmFsdWF0aW9uIGZvciBzaW5nbGUgc3RlcA0KYGBge3J9DQpsaWtlbGlob29kX3N0ZXAgPC0gZnVuY3Rpb24odGhldGEsIHgsIHhfbmV4dCwgZGVsdGF0KXsNCiAgI0RlZmluaW5nIHR1bmFibGUgYW5kIGNvbnN0YW50IHBhcmFtZXRlcnMNCiAgbGFtYmRhIDwtIHRoZXRhWzFdDQogIE4xIDwtIHRoZXRhWzJdDQogIEMgPC0gdGhldGFbM10NCiAgc2lnbWExIDwtIHRoZXRhWzRdDQogIHNpZ21hMiA8LSB0aGV0YVs1XQ0KICBzaWdtYTMgPC0gdGhldGFbNl0NCiAgDQogICNEZWZpbmluZyBkcmlmdCBmdW5jdGlvbg0KICBkcmlmdCA8LSBmdW5jdGlvbih4KXsNCiAgICB4MSA8LSB4W1sxXV0NCiAgICB4MiA8LSB4W1syXV0NCiAgICB4MyA8LSB4W1szXV0NCiAgICB5IDwtIG51bWVyaWMobGVuZ3RoID0gMykNCiAgICB5WzFdIDwtIGxhbWJkYSAtIHJobyp4MSAtIGV0YSp4MSp4Mw0KICAgIHlbMl0gPC0gZXRhKngxKngzIC0gZGVsdGEqeDINCiAgICB5WzNdIDwtIE4xKmRlbHRhKngyIC0gQyp4Mw0KICAgIHkNCiAgfQ0KICAjIG5lZWQgdG8gYXNzaWduIHN1YnBhcmFtZXRlcnMgb2YgdGhldGEgdG8gdGhlIG5hbWVzIHVzZWQgaW4gdGhlIGNvZGUNCiAgeiA8LSBudW1lcmljKDMpDQogIHpbMV0gPC0geFtbMV1dDQogIHpbMl0gPC0geFtbMl1dDQogIHpbM10gPC0geFtbM11dDQogICN6X25leHQgPC0gbnVtZXJpYygzKQ0KICAjel9uZXh0WzFdIDwtIHhfbmV4dA0KICAjYnJvd3NlcigpDQogIHkgPC0gZHJpZnQoeCkNCiAgYSA8LSBsb2coKCgyKnBpKmRlbHRhdCleMykgKiAoelsxXSp6WzJdKnpbM10qc2lnbWExKnNpZ21hMipzaWdtYTMpXjIpDQogIGIgPC0gKCh4X25leHQgLSB4ICsgeSpkZWx0YXQpXjIpKmMoKHpbMV0qc2lnbWExKV4oLTIpLCAoelsyXSpzaWdtYTIpXigtMiksICh6WzNdKnNpZ21hMyleKC0yKSkvZGVsdGF0DQogIGMgPC0gYltbMV1dKmJbWzJdXSpiW1szXV0NCiAgYStjDQogIH0NCmBgYA0KDQpgYGB7cn0NCiNUZXN0IGNodW5rDQp0aGV0YSA8LSBjKGxhbWJkYSwgTjEsIEMsIHNpZ21hMSwgc2lnbWEyLCBzaWdtYTMpDQpsaWtlbGlob29kX3N0ZXAodGhldGEgPSB0aGV0YSwgeCA9IHhfY29hcnNlWzEsXSwgeF9uZXh0ID0geF9jb2Fyc2VbMixdLCBkZWx0YXQgPSAwLjEpDQojV2UgZ2V0IGEgbnVtYmVyLCBzbyB0aGF0J3MgZ29vZC4NCmBgYA0KSW1wbGVtZW50aW5nIGxpa2VsaWhvb2QgZm9yIGVudGlyZSBwcm9jZXNzDQpgYGB7cn0NCmxpa2VsaWhvb2RfdG90YWwgPC0gZnVuY3Rpb24odGhldGEsIHByb2Nlc3MsIGRlbHRhdCl7DQogIG5vc3RlcHMgPC0gbGVuZ3RoKHByb2Nlc3NbLDFdKS0xDQogIGN1bXVsYXRpdmVfbGlrZWxpaG9vZCA8LSAwDQogIGZvciAoaSBpbiAoMTpub3N0ZXBzKSl7DQogICAgY3VtdWxhdGl2ZV9saWtlbGlob29kIDwtIGN1bXVsYXRpdmVfbGlrZWxpaG9vZCArIGxpa2VsaWhvb2Rfc3RlcCh0aGV0YSA9IHRoZXRhLCB4ID0gcHJvY2Vzc1tpLF0sIHhfbmV4dCA9IHByb2Nlc3NbaSsxLF0sIGRlbHRhdCA9IGRlbHRhdCkNCiAgfQ0KICBjdW11bGF0aXZlX2xpa2VsaWhvb2QNCn0NCmBgYA0KVGVzdGluZyB0b3RhbCBsaWtlbGlob29kIGZ1bmN0aW9uDQpgYGB7cn0NCmxpa2VsaWhvb2RfdG90YWwodGhldGEgPSB0aGV0YSwgcHJvY2VzcyA9IHhfY29hcnNlLCBkZWx0YXQgPSAwLjEpDQojV2UgZ2V0IGEgbnVtYmVyLCB3aGljaCBpcyBwcm9taXNpbmcNCmxpa2VsaWhvb2RfdG90YWwodGhldGEgPSB0aGV0YV8wLCBwcm9jZXNzID0geF9jb2Fyc2UsIGRlbHRhdCA9IDAuMSkNCmBgYA0KV2Ugc2VlIHRoYXQgZXZhbHVhdGluZyB0aGUgdG90YWwgbGlrZWxpaG9vZCBmdW5jdGlvbiByZXR1cm5zIHRoZSBzYW1lIG51bWJlciBmb3IgdHdvIGRpZmZlcmVudCBzZXQgb2YgcGFyYW1ldGVycy4NClRoaXMgaW5kaWNhdGVzIHRoYXQgdGhhdCB0aGUgdG90YWwgbGlrZWxpaG9vZCBmdWNudGlvbiBkb2Vzbid0IGluIGZhY3QgdXNlIHRoZSBwYXJhbWV0ZXIgcHJvdmlkZWQgZm9yIGV2YWx1YXRpb24uDQpJdCBzZWVtcyB0aGF0IGxpa2VsaWhvb2Rfc3RlcCB0YWtlcyBpdHMgcGFyYW1ldGVycyBmcm9tIHRoZSBnbG9iYWwgZW52aWVyb25tZW50IHJhdGhlciB0aGFuIGl0cyBhcmd1ZW1lbnQuDQpUaHMgd291bGQgZXhwbGFpbiB3aHkgdGhlIGxpa2VsaWhvb2QgaXMgZXZhbHVhdGVkIHRvIHRoZSBzYW1lIGZvciBhbnkgcGFyYW1ldGVyLg0KTm93IHdlIGdvdCBkaWZmZXJlbnQgbnVtYmVycy4NCg0KTm93LCB3ZSB0cnkgdG8gb3B0aW1pemUgdGhpcyBmb3Igc29tZSBwYXJhbWV0ZXIgdmVjdG9yIHRoZXRhLg0KYGBge3J9DQp0aGV0YV8wIDwtIGMoNTAsIDIwMDAsIDEwLCAwLjUsIDAuNSwgMC41KQ0Kb3B0aW0ocGFyID0gdGhldGFfMCwgZm4gPSBsaWtlbGlob29kX3RvdGFsLCBwcm9jZXNzID0geF9jb2Fyc2UsIGRlbHRhdCA9IDAuMSwgbG93ZXIgPSBjKDAsMCwwLDAsMCwwKSwgbWV0aG9kID0gIkwtQkZHUy1CIikNCmBgYA0KSXQgd291bGQgYXBwZWFyIHRoYXQgaSBoYXZlIG1hZGUgYSBtaXN0YWtlOiBPcHRpbSBvdXRwdXRzIHRoZSBpbml0aWFsIHBhcmFtZXRlcnMsIGFuZCBzdGF0ZXMgdGhhdCBjb252ZXJnZW5jZSB3YXMgbm90IHJlYWNoZWQuDQpBbm90aGVyIG1pc3Rha2UgaXMgbm90IGZpbmlzaGluZyBhIGNvZGluZyBwcm9qZWN0IGJlZm9yZSB0YWtpbmcgYSAyLXdlZWsgYnJlYWs6IEl0J3MgYSBsaXR0bGUgZGlmZmljdWx0IHRvIGp1bW8gaW4gYWdhaW4uDQoNCkl0IHRha2VzIHNvbWUgdGltZSB0byBydW4gbm93LCBhbmQgd2UgZG9uJ3QgZ2V0IGNvbnZlcmdlbmNlIG9yIHRoZSByaWdodCByZXN1bHRzLg0KVGhlIHRvdGFsIGxpa2VsaWhvb2QgZnVuY3Rpb24gd2FzIGNhbGxlZCBvdmVyIDQwMCB0aW1lcywgcHJvYmFibHkgYmVjYXVzZSBvZiBsYWNrIG9mIGNvbnZlcmdlbmNlLCB3aGljaCBtaWdodCBhbHNvIGJlIHdoeSBpdCB0YWtlcyBzbyBsb25nLg0KUXVlc3Rpb246IGFyZSB3ZSBtYXhpbWl6aW5nIG9yIG1pbmltaXppbmc/DQpXaGF0IGlzIHRoZSBwYXJhbWV0ZXIgc3BhY2U/DQpXaGF0IGV4YWN0bHkgaXMgdGhlIHF1YW50aXR5IGJlaW5nIGNvbXB1dGVkIGluIHRoZSBzdGVwIGxpa2VsaWhvb2QgZnVuY3Rpb24/IEkgdGhpbmsgaXQgaXMgdGhlIGxpa2VsaWhvb2QsIGZvcm11bGEgZnJvbSB0aGUgc2xpZGVzLg0KTGV0J3MgdHJ5IHdpdGggc29tZSBwYXJhbWV0ZXJzIHRoYXQgYXJlIGNsb3NlciB0byB0aGUgcmVhbCBvbmVzLg0KV2UgZ2V0IGNvbXBsZXRlbHkgYXRyb2Npb3VzIGVzdGltYXRlcy4NCldlIHRyeSBib3VuZGluZyBiYXJhbWV0ZXIgdmFsdWVzIGZyb20gYmVsb3cgd2l0aCB6ZXJvIGFuZCB1c2luZyBMLUJGR1MtQiBtZXRob2QuDQpJdCBzZWVtcyB0aGF0IHdlIGdldCBpbmZpbml0ZSB2YWx1ZXMgZnJvbSB0b3RhbCBsaWtlbGlob29kLCB3aGljaCBpcyB1bmZvcnR1bmF0ZS4gTWF5YmUgdGhlIGxvZyBpcyBmdWNraW5nIG1lPw0KDQoNCg0K